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Abstract 

We study the impact of competing time delays in coupled stochastic syn- 
chronization and coordination problems. We consider two types of delays: 
transmission delays between interacting elements and processing, cognitive, 
or execution delays at each element. We establish the scaling theory for 
the phase boundary of synchronization and for the steady-state fluctuations 
in the synchronizable regime. Further, we provide the asymptotic behavior 
near the boundary of the synchronizable regime. Our results also imply the 
potential for optimization and trade-offs in synchronization problems with 
time delays. 
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1. Introduction 

Coordinating, distributing, and balancing resources in networks is a com- 
plex task as these operations are very sensitive to time delays (J, 0]. To 
understand and manage the collective response in these coupled interacting 
systems, one must understand the interplay of stochastic effects, network 
connections, and time delays. In synchronization, coordination, and con- 
sensus problems in coupled interacting systems [l|-l5|, individual units at- 
tempt to adjust their local state variables (e.g., pace, load, orientation) in 
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a decentralized fashion. They interact or communicate only with their local 
neighbors in the network, often with explicit or implicit intention to improve 
global performance. Applications of the corresponding models range from 
physics, biology, computer science to control theory, including synchroniza- 
tion problems in distributed computing 0, 0, 0], coordination and control 
in communication networks [3, 0, 0-12], flocking animals 1315], bursting 



neurons 



16r|l9(, and cooperative control of vehicle formation [20| 



In this Letter, we study the impact of competing, finite non-zero time 



delays in stochastic synchronization or coordination prob 



present in most real communication, information (3,0,0, 10 1 



logical systems [24M26I] including neurobiological networks [17H19 



cms 
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which are 
23| , and bio- 
(Through- 



out this Letter, we use the terms coordination and synchronization synony- 
mously.) Delays can be attributed to both non-zero transmission times be- 
tween the nodes and to non-zero finite times it takes to process (possibly in- 
cluding cognitive delays) and execute the desired action at the nodes. Here, 
we investigate the importance and impact of these two types of delays in 
a simple synchronization problem in noisy environment with two linearly 
coupled nodes. 

Singularities in critical phenomena and phase transitions [27|, which are 
often present in coupled interacting systems consisting of a large number 
of nodes N, are typically associated with progressively more eigenvalues of 
the coupling operator (e.g., the Laplacian) getting arbitrarily close to zero. 
Strictly speaking, these singularities are exhibited only by systems approach- 
ing the thermodynamic limit (where the density of eigenvalues does not van- 
ish sufficiently fast, or itself becomes singular about zero in the iV— s-oo limit). 
For example, in spatially-embedded physical systems these singularities are 
typically exhibited by the relevant response functions and fluctuations in 
the long- wavelength limit [4]. In complex networks 28-31] these singulari- 
ties can be suppressed as a result of sufficient amount of randomness in the 



connectivity pattern [4j, I32N35 



In contrast, the instability governed by time delays is associated with a 
single mode exceeding a threshold value (in a simple case, associated with 
the eigenmode of the network Laplacian with the largest eigenvalue 0, fiol] ) ■ 
Therefore, the underlying instability is present even in the simplest network 
with two nodes (N=2). Here we focus on a two- node network, which qual- 
itatively captures the generic features of the coordination behavior when 
the delays are present due to both transmission between nodes and process- 
ing/execution at each node. For simplicity, we will refer to this instability 
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as "critical", even though this singularity does not require infinitely many 
degrees of freedom. In networks, consisting of a large number of nodes, the 
effect of time delays is not qualitatively different, but can be "amplified" by 
heterogeneous 



(or scale- 
of uniform time delays [l|. 




connectivity patterns: in the case 



the effective coupling of the most relevant 
singular mode is the largest eigenvalue of the coupling operator, which itself 
can diverge with the system size (3, 36 38], severely limiting synchronizability 
and coordination. 



2. A Stochastic Model with Local and Transmission Delays 



Differential equations with delays [39] describing complex systems have 
a long history, originally motivated by the emergence of business and eco- 



nomics cycles [4CH42| . and also naturally appearing in the context of stability 



of ecologic al sy stems, in models in population dynamics, and in game theory 



24j. l25l. |43|-|46|. There have been recent works combining stochastic differen- 
tial equations with delays 47H49| with applications ranging from population 
dynamics, epidemiology, and immunology to cell kinetics and finance 50|, |51 



Here we consider a model for local coordination where time delays are 
attributed to two separate origins: one is the transmission between the two 
nodes, the other is processing the information and executing the action at 
each node, denoted by r tr and r Q , respectively. We investigate the simplest 
stochastic model where the coordination or synchronization attempt between 
the two nodes, in terms of the relevant state variables hi, is captured by linear 
relaxation 



d t hx(t) = -A[/li(t-T )-/l 2 (t-To-Tt r )]+77i(t) 

dMt) = -A[fc 2 (t-T )-Mt-T -T tr )] + r72(t) • (!) 

Here, rji is delta correlated noise with (77$) = and (Vi{t)Vj{t')) = 2D5ij8(t—1f) 
with noise intensity D, i,j — 1,2. A > is the coupling strength between 
the two nodes. For initial conditions, we use hi(t)=0 for t<0. 

To simplify notation we introduce r = r + r tr and 7 = r Q / (r + r tr ) = t /t 
(0 < 7 < 1). Further, since we are interested in the synchronization (or 
coordination) between the two nodes, we focus on the relative difference 
u{t) = h 2 (t) — hi(t) which is governed by 

d t u{t) = -Xu(t - 7T) - Xu(t - t) + , (2) 
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where (£) = and (£(£)£(£')) = 4D5(t— £'). The special case 7=1 of the above 
equation has been investigated in our earlier work [2[. Here, we study the im- 
pact of both types of delays, corresponding to the general case 0<7<1. Our 
quantity of interest is (u 2 (t)), capturing the relative deviation of the relevant 
state variables on the two nodes. By definition, the system is synchronizable 
if the fluctuations reach a finite steady state, (u 2 (oo)) < 00. In the absence 
of time delays (t=0) one immediately finds (u 2 (t)) = (D/X)(l — e~ iXt ) (52| . 
i.e., the system is synchronizable for any A > 0. Further, the stronger the 
coupling, the better the synchronization: (w 2 (oo)) = D/X is a monotonically 
decreasing function of A. 

Next we study and analyze the case with time delays. Employing stan- 
dard Laplace transform (2, 53], one can immediately write the formal solution 
for Eq. © 

where s a , a = 1, 2, . . ., are the zeros of the characteristic equation 

g(s) = s + Ae" 7rs + Xe~ TS = (4) 
on the complex plane. Then for the noise- averaged fluctuations one finds 

= E v ? v 1 \ 



-4D(1 - e (sQ+s ' 3)t ) 



7Are"T rs « - Are- rs ")(l - 7Are~ 7TSs - XTe~ TS P){s a + sp) 



—4Dt(1 — e ( ~ Za+z ^ T ) 



. _ - 7Ae"T 2a - Ae" 2 «)(l - -fAe" 1 ^ - ke~ z P)(z a + zp) 



(5) 



where in the last expression of the above equation we introduced the scaled 
variables z a = rs a and A = Ar. From Eq. (T4]) and from the definition of 
these scaled variables it is evident that the solutions of the scaled 

characteristic equation 

z + Ae" 7 * + Ae z = , (6) 

and consequently, the solutions depend only on A, (A). From 

the structure of the above characteristic equation it follows that if z is a 
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Figure 1: (a) Time series (u 2 (t)) for r=1.00, 7=0.50 for different values of the coupling 
constant A. Here, and throughout this paper, D=\ and Ai=0.01. (b) Synchronizability 
threshold in terms of the scaled variable At vs 7. Data points were obtained by numerically 
integrating Eq. ([2]) [54 ] . The solid line represents the exact analytic expression Eq. ([7]). 



solution of Eq. so is its complex conjugate z*. From Eq. ([5]) it is clear 
that synchronization can only be achieved if Re(z Q )<0 for all a. To identify 
the boundary of the region of synchronizability, one has to find the solu- 



tion^) with a vanishing real part, i.e., z = x + iy with x=0 25|, |40|, |42|, |44 
Elementary analysis yields yf=±n/ (1 + 7) and 

Thus, for a fixed 7, the system is synchronizable if < Ar < A c (7). Results 
obtained by numerically integrating Eq. (j2J) together with the analytic 
expression Eq. (|7|) are shown in Fig. [TJ We will discuss the phase diagram 
and some limiting cases in terms of the original variables, the local delay r Q 
and the transmission delay r tr , in the final section. 



3. Scaling and Asymptotics in the Steady State 

Now we turn to analyzing the steady-state fluctuations, in particular, 
their scaling behavior in the synchronizable regime, < A < Ac (7). Here, 
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the fluctuations remain finite, and in the steady state (t— >oo) from Eq. (151) 
one obtains 

(u 2 (oo)} = Drf( 1 ,A) , (8) 

where 

f(ry A) = — (9) 

w ' y (1 -7Ae-T*» - Ae" 2 «)(l -^Ae~^ - Ae- Z ?)(z a + Zp) w 

is the scaling function for the steady-state fluctuations. [Recall that z a = 
z a (A) are the solutions of the scaled characteristic equation Eq. (Q.} Thus 
for a given 7, 

^ = /(At) , (10) 

where in our notation, we suppressed the 7 dependence to highlight the scal- 
ing behavior of the fluctuations, valid for each 7 separately. Figure [2] shows 
the steady-state fluctuations before (a) and after (b) scaling, and demon- 
strates the data collapse for the scaled variables according to Eq. (fTUj) . The 
scaling function /(A) is a non-monotonic function of its argument, diverging 
at A=0 and A=A C (7), and exhibiting a single minimum between these points 
[Fig. [2(b)]. This non-monotonic feature of the scaling function with a single 
minimum between < At < A c (7) is present for all < 7 <1 [Fig. [3]. Thus, 
for fixed non-vanishing and finite delays, there is an optimal value of the 
coupling constant A for which the steady-state fluctuation attains its mini- 
mum value. For stronger couplings, the overall coordination between the two 
nodes weakens, and for A > A c (j)/t, it completely deteriorates. 

Next, we briefly discuss the asymptotic behavior of the scaling func- 
tion near the boundaries of the synchronizable regime. The fluctuations 
of (u 2 (oo)) diverge at the end points of this interval [as at least for one a, 
Ke(z a ) 0], indicating the breakdown of synchronization. Near these end- 



points, the sum in Eq. (Q is dominated by the term(s) where Re(z a ) ~ [55 
These are the solutions which have (negative) real parts with the smallest 
amplitude. As we show in Appendix A, to leading order, 



as A— ^0, and 



/(A) =4 <") 
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Figure 2: (a) Steady-state fluctuations as a function of the coupling strength A for the 
various delays for 7=0.5. Data points are obtained by numerically integrating Eq. ([2]). 
(b) Same data as in (a) scaled according to Eq. (jTU)) . {u 2 (oo))/r vs At. The dashed lines 
represent the asymptotic behaviors of the scaling function near the two endpoints of the 
synchronizable regime [H3|, Eqs. (fTTj) and (fl2j) . respectively, while the solid line (running 
precisely through the data points) represents the full approximate scaling function /(At), 
Eq. (13). 



as A— >A C (A<A C ) with 01(7) given in Appendix A [Eq. (1A.7j) ]. From the 



numerical results [Fig. EJ^b)] it is also apparent that the scaling function 
varies slowly between (and away from) the singular points, thus, /(A) can be 
reasonably well approximated [55| throughout the full synchronizable regime 
0<A<A c ( 7 )by 

with c 2 (7) also given in Appendix A [Eq. (lA.llj) . 

Figure EJ^b) and Fig. [3] show that the above approximate scaling function 
Eq. (|T3|) (being asymptotically exact near the singular points) matches the 
numerical data very well. In particular, it captures the basic non-monotonic 
feature of the results obtained from numerical integration, exhibiting a single 
minimum 

A mi „(7) = —^7= (14) 
in the < A < A c (7) interval. 
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Figure 3: Scaled steady-state fluctuations (u 2 (oo))/r vs At for various 7 values. Data 
points are obtained by numerically integrating Eq. @. Solid lines represent the full 
approximate scaling function /(Ar) for each 7, Eq. fj 13[) . 



As can also be seen in Fig. [31 the theoretical asymptotic behavior, cap- 
tured by the approximate scaling function Eq. f[T5j) becomes less accurate for 
small 7 near A c (7). Other than lacking higher-order corrections to the asymp- 
totic expressions, this is due in part to the time discretization in the numerical 
integration |54|. For sufficiently small 7 values, the condition At « 7T will 



not hold, and deviations between the results of the time-discretized numer- 
ical scheme and those of the continuous system Eq. ([2]) will become more 
significant and noticeable. 

4. Discussion and Summary 

Having established the scaling theory for the phase boundary [Eq. (j7])] 
and for the fluctuations [Eq. ffl3]) ]. it is insightful to express our main find- 
ings explicitly in terms of the two types of delays appearing in the original 
formulation of the problem, Eq. ([T]). From Eq. ([7]), for the boundary of the 
synchronizable regime one immediately finds 

A(2r D + 7tr) = r . (15) 
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Figure 4: Synchronizability phase diagram on the Art r -AT plane [Eq. (Tl"5j) ]. The shaded 
area indicates the synchronizable regime. The boundary of this region approaches the 
horizontal line Ar =l/2 in the limit of Art r — too. Further, At =7t/4 when Art r =0. 

While explicitly expressing the critical line r Q vs r tr is prohibitive due to the 
implicit nature of Eq. (Tl5|) . one can produce a plot for it numerically without 
any difficulties [Fig. H] . 

We can gain further insight into the different impact of the two types 
of delays by considering two limiting cases. First, consider the case when 
r o/ r tr << 1; i-e., when the transmission delays are much larger than the 
local processing, cognitive, or execution delays. This is equivalent to the 
7<<1 limit in our scaling expressions. From Eq. (J7|) one finds A c ~l/27 or 
(At ) c =1/2. Thus, there is no singularity in the fluctuations for any finite r tr 
provided that Ar G < 1/2. Further, from Eq. ( fl~3l) (with the coefficients given 
in Appendix A) for the steady-state fluctuations in the same limit we find 

In the other limiting case, r tr /r << 1, i.e., the transmission delays are 
much smaller than the local processing delays. This is equivalent to the 7—7-1 
limit in our scaling expressions. In this limit Eq. (J7|) reduces to A c ~7r/4 or 
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(At q ) c = 7r/4. The steady-state fluctuations approach 



A [ 7r 7r/4 — Ar D 



7T 



provided that At q < 7r/4. 

Figure HI and Eqs. (fT6j) and (ITTj) highlight the subtle differences between 
the impacts of the two types of delays. We may regard the local delays r D as 
the dominant ones, in that as long as Ar D < 1/2, there are no singularities for 
any finite r tr , and (u 2 (oo)) increases linearly with r tr as r tr — t-oo [Eq. (JIB])]. 
On the other hand, for every r tr , there is a sufficiently large r Q such that the 
fluctuations become singular. In particular, when the transmission delays 
are much smaller than the local processing delays, the fluctuations diverge 
as At ^tt/4 [Eq. ([T7j)]. 

Inside the synchronizable regime, for fixed r D and r tr (with the exception 
of r o =0 [56]), the steady-state fluctuations (w 2 (oo)) always exhibit a single lo- 
cal minimum as a function of the coupling constant A [Eq. Q14j) ] . This feature 
naturally presents scenarios for optimization and trade-offs. In particular, in 
real systems, the effective coupling constant between two interacting nodes 
can be controlled by the frequency of pairwise communications. This im- 
plies that too much communication can cause more harm than good, and 
further, there is an optimal rate of communications between the nodes which 
minimizes the size of the fluctuations. Also, as was already shown for the 
special 7=1 case |2|, in large network-coupled systems, decreasing connectiv- 
ity may be beneficial if the system is too close or beyond the synchronization 
boundary. 

In this Letter we only considered linear couplings between two nodes 
in the presence of noise and competing time delays. We established the 
phase diagram and the scaling theory for synchronizability, and provided 
the asymptotic behavior for the relevant scaling functions. Nonlinearities 
are undoubtedly important in real systems, and will likely further increase 
the complexity of the already rich phase diagram, such as recurring and 



alternating patterns of synchronizable and unsynchronizable regions [18|, [19 



22 . 
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Appendix A. Asymptotic Behavior of the Scaling Function Near 
the Synchronization Thresholds 



Here, we employ the method in Ref. [55J to calculate the dominant con- 
tributions in Eq. ([9]) near the boundaries of the synchronizable regime. We 
assume that solutions of the characteristic equation 

z + Ae" 7z + Ae~ z = (A.l) 

change continuously with the parameter A. Thus, if z=z Q is a solution for 
A=A G , then for a small change in the parameter, A = A + SA, the corre- 
sponding solution can be written as z = z Q + 5z. Substituting this into the 
characteristic equation, to lowest order we find 

p T^o _|_ p Zo 

Sz ~ - - SA + O ((SA) 2 ) . (A.2) 

For A=0, there is a single solution with vanishing real part, z—0, thus for 
small A 

z(A) ~ -2A + 0(A 2 ) . (A.3) 

The dominant contribution for the scaling function as A-+0 comes from the 
corresponding term in Eq. ([9]), to leading order yielding 

/(A) ~ ~ 4 AX = 4 . (A.4) 
JK ' 2(-2A) A v ' 

For A=A C (7) [Eq. there is a pair of solutions (complex conjugates) 
with vanishing real parts z = ±i^-. When A is in the vicinity of A c (A ~ 
A c + SA), to lowest order, these solutions behave as 

e =R7?/c _|_ e T*2/c 

z±(A)~±iy c -- : — 5 A , (A. 5) 
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where y c = . The dominant contributions for the scaling function as 
A— i-A c then come from the two terms in Eq. ([9]) when (a = ±, (3 = =)=), 
yielding 

-8 



/(A) 



1 - 7A c e-^ c - A c e"^)(l - 7A c e^ c - A c e^)(^+ + zJ) 
ci(7) 



A c (7) - A 



(A.6) 



where 



Ci(t) = 1 • (A.7) 

(1 + 7)A C + (1 + 7)A C cos(vr^) - cos(^) - cos(^) 

Finally, we obtain the approximate scaling function for the full 0<A<A C (7) 



interval, using some heuristics following Ref. 55]. As observed from our nu- 
merical results [Fig. [2] , the scaling function varies slowly between (and away 
from) the two singular points. Then, it can be approximated by 

In principle, the constant 02(7) could be determined by matching the mini- 
mum value of the scaling function. Since it is not known analytically, instead 



we resort to the heuristics of Ref. [55] where the constant 02(7) is determined 
in such a way that it matches next-to-leading order corrections of the asymp- 
totic behavior, e.g., near A=0. To that end, we find the next-to-lowest order 
corrections to the solution of Eq. (lA.ljl in the vicinity of A=0, 

z{A) ~ -2A- 2(1 + 7 )A 2 + 0(A 3 ) . (A.9) 

Keeping the relevant orders in the dominant term in Eq. (Q, we obtain 

-4 

f(A) ~ 

JV } (1- 7 A-A) 2 2(-2A-2(1 + 7 )A 2 ) 

1 

~ [1-(1+7)A] 2 A(1 + (1+ 7 )A) 

~ i + (l + 7 ). (A.10) 

In order to match this next-to-leading order correction as A— >0 with the 
proposed approximate scaling function Eq. flA.8j) . one must have 



= (A - U) 
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